clear all
tic
n=20;

%% Constrained states

Adj=1*diag(ones(n-1,1),1)+1*diag(ones(n-1,1),-1);Adj(end,1)=1;Adj(1,end)=1;
a=[0,1];
for k=2:n
    free=a;
    for j=1:k-1;
        if Adj(k,j)==1
            free=free(bitget(free,j)==0);
        end
    end
    a=[a,free+2^(k-1)];
end
toc
 
%% sectors with differnt symmetry

% following http://physics.bu.edu/~sandvik/vietri/dia.pdf

tic
[R,M,s,ind]=checkstate_1D(a,n);
toc
%% Hamiltonain construction

%% X-term 
tic
Sx=cell(n,1); % Sx is the sum of local x pauli matrices
for k=1:n
    Sx{k}=sparse(M{k},M{k});
end
for j=1:n
    
    sout=bitxor(s{n},2^(j-1)); %%?????? why s{n}????? Because it contains all the (contrained allowed) representatives
    [repout,lout]=representative(sout,n);
    
    for k = 1:n
    r=repout(ind{k});
    l=lout(ind{k});
    [Q,S]=ismember(r,s{k});
    
    index1=find(S>=1);%mapped from
    index2=S(index1);%mapped to
    
    Ra=R{k}(index1);Ra=Ra(:);
    Rb=R{k}(index2);Rb=Rb(:);
    lb=l(index1);lb=lb(:);
    elems=sqrt(Ra./Rb).*exp(1i*2*pi*k*(lb)/n);
    
    Sx{k}=Sx{k}+sparse(index2,index1,elems,M{k},M{k});
   
    end
end

for k=1:n 
    Sx{k}=(Sx{k}+Sx{k}')/2;
end
toc
%% Z-term 
tic
Pe=cell(n,1); % Sx is the sum of local x pauli matrices
for k=1:n
    ne=zeros(size(s{k},1));
    for j=1:n
        ne=ne+bitget(s{k},j);
    end
    Pe{k}=sparse(1:M{k},1:M{k},ne);
end
toc

%% Correlations
PePe=cell(n,n-1); % This is Sum(i) Pe(i)Pe(i+d)
for k=1:n
    for d=1:n-1
        nene=zeros(size(s{k},1));
        for j=1:n
            nene=nene+bitget(s{k},j).*bitget(s{k},mod(j+d-1,n)+1);
        end
        PePe{k,d}=sparse(1:M{k},1:M{k},nene);
    end
end

%%%%%
psi0=eye(size(s{n},1));
tic
[PSI,tvals]=expmv_tspan(-1i*Sx{n},psi0,0,10,100);
toc